Code
library(dplyr)
library(tidyr)
library(vegan)
library(ggplot2)
library(targets)
library(ggrepel)
library(patchwork)
theme_set(theme_light())Appendix S1
This document contains all data analyses presented in the main text as well as supporting information mentioned and extra information on the analysis for those interested (not cross-referenced).
The analysis was organized using the make-line pipeline tool provided by targets (Landau 2021). You will find raw data as well as all functions used in the pipeline in the project’s research compendium created with the rcompendium package (Casajus N. 2022).
You can take a look at the pipeline in order to have an idea of the different steps for creating the main models and figures. Each target is defined at the _targets.R file in the research compendium.
Analysis pipeline. Total runtime and storage size of each object (target) are shown.
| Supporting | Main text |
|---|---|
| Figure S 1 | Figure 1 |
| Figure S 2 | Figure 2 |
| Figure S 3 | Figure 3 |
| Figure S 4 | Figure 4 |
| Figure S 5 | Figure 5 |
| Figure S 6 | Figure 6 |
| Figure S 7 | Figure 7 |
| Figure S 9 | Figure 8 |
| Figure S 8 |
Unfortunately links only redirect to tabsets open by default but not hidden ones. Please look for tables related to analyses by phylum or by compartment in their respective tabsets.
| Supporting Tables |
|---|
| Table S 1 |
| Table S 2 |
| Table S 3 |
| Table S 4 |
| Table S 6 |
| Table S 7 |
| Table S 8 |
| Table S 9 |
| Table S 10 |
| Table S 11 |
| Table S 12 |
| Table S 13 |
| Table S 14 |
library(dplyr)
library(tidyr)
library(vegan)
library(ggplot2)
library(targets)
library(ggrepel)
library(patchwork)
theme_set(theme_light())Take a look at the structure of the final complexity dataset
tar_load(comp)
knitr::kable(head(comp))
| Mini | Site | Point | Sample | Branching_density | Broken_density | D_bin | D_gray | L | I | S | Total_Density | Dry_weight | Sphericity | DR1 | DR2 | DR3 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MBE_1 | Belle-ile | Belle-ile 1 | 1 | 5.8 | 0.0 | 1.1231 | 2.6562 | 12.85333 | 9.143333 | 8.303333 | 181696 | 183.61 | 0.8371338 | 0.6460062 | 0.8153846 | 0.7113589 |
| MBE_1 | Belle-ile | Belle-ile 1 | 2 | 1.0 | 0.0 | 1.0654 | 2.6463 | 16.07667 | 7.883333 | 5.393333 | 181696 | 183.61 | 0.6122605 | 0.3354758 | 0.7669267 | 0.4903587 |
| MBE_1 | Belle-ile | Belle-ile 1 | 3 | 3.0 | 0.2 | 1.0722 | 2.6450 | 15.30000 | 8.103333 | 7.190000 | 181696 | 183.61 | 0.7470808 | 0.4699346 | 0.8873818 | 0.5296296 |
| MBE_1 | Belle-ile | Belle-ile 1 | 4 | 4.8 | 0.4 | 1.1078 | 2.6492 | 12.91667 | 8.476667 | 6.203333 | 181696 | 183.61 | 0.7057078 | 0.4802581 | 0.6613704 | 0.6562581 |
| MBE_1 | Belle-ile | Belle-ile 1 | 5 | 4.2 | 0.0 | 1.1305 | 2.7019 | 24.54667 | 14.523333 | 8.970000 | 181696 | 183.61 | 0.6088477 | 0.3654264 | 0.6434838 | 0.5916621 |
| MBE_1 | Belle-ile | Belle-ile 1 | 6 | 3.0 | 0.4 | 1.1198 | 2.6641 | 16.35667 | 9.416667 | 6.783333 | 181696 | 183.61 | 0.6684949 | 0.4147137 | 0.7249304 | 0.5757082 |
Most variables are not normally distributed so we apply a box-cox (Box and Cox 1964) transformation and check correlations between the transformed variables
tar_read(pairscomp)Let’s check it without removing highly correlated variables
tar_load(pcatotal)
tar_read(sptot)We had up to 6 relevant axes to look at, but we know we have several highly correlated variables
tar_read(pcatot)tar_read(pcatotsel)Call: rda(X = num_sel, scale = TRUE)
Inertia Rank
Total 5
Unconstrained 5 5
Inertia is correlations
Eigenvalues for unconstrained axes:
PC1 PC2 PC3 PC4 PC5
1.6866 1.1941 0.8642 0.6783 0.5768
tar_read(sptotsel)Let’s check the biplots of the selected variables
tar_read(pcasel)tar_read(spmedsel)Let’s check the biplots of the selected variables
Figure 1
tar_read(pcamed)Let’s check if we find similar groupings using a ward’s cluster
tar_read(clustcomp)tar_load(corcomp)
plot(corcomp[[1]])
plot(corcomp[[2]])
tar_read(cortestcomp1)
tar_read(cortestcomp2)
Pearson's product-moment correlation
data: med_sel2$PC1_c and med_sel2$PC1_score
t = -16.404, df = 28, p-value = 6.824e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.9769962 -0.9000432
sample estimates:
cor
-0.9517122
Pearson's product-moment correlation
data: med_sel2$PC2_c and med_sel2$PC2_score
t = 25.364, df = 28, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.9557115 0.9900324
sample estimates:
cor
0.978924
tar_load(faunadata)
knitr::kable(head(faunadata))
| Point | Site | Replicate | Method | Habitat | Year | Season | Date | Species | Abundance |
|---|---|---|---|---|---|---|---|---|---|
| Belle-ile 1 | Belle-ile | 1 | Grab | Maerl | 2007 | Spring | 2007-03-11 00:00:00 | Abludomelita gladiosa | 1 |
| Belle-ile 1 | Belle-ile | 1 | Grab | Maerl | 2007 | Spring | 2007-03-11 00:00:00 | Abra alba | 1 |
| Belle-ile 1 | Belle-ile | 1 | Grab | Maerl | 2007 | Spring | 2007-03-11 00:00:00 | Ampharete spp. | 1 |
| Belle-ile 1 | Belle-ile | 1 | Grab | Maerl | 2007 | Spring | 2007-03-11 00:00:00 | Amphipholis squamata | 2 |
| Belle-ile 1 | Belle-ile | 1 | Grab | Maerl | 2007 | Spring | 2007-03-11 00:00:00 | Anapagurus hyndmanni | 1 |
| Belle-ile 1 | Belle-ile | 1 | Grab | Maerl | 2007 | Spring | 2007-03-11 00:00:00 | Aonides oxycephala | 3 |
samplespoint <- faunadata %>%
mutate(Date=as.factor(format.Date(faunadata$Date,"%Y"))) %>%
distinct(Site, Year, Point, Replicate) %>%
group_by(Site, Year, Point) %>%
mutate(Year = as.factor(Year)) %>%
summarise(count = n()) %>%
ungroup() %>%
mutate(State = factor(case_when(count==3 ~ "Good", count < 3 ~ "Sample(s) missing", count > 3 ~ "Too many samples")))
ggplot(samplespoint, aes(x=Year,y=Point,fill=State)) +
geom_tile(alpha=0.5) +
scale_fill_manual(values=c("palegreen3","lightblue", "indianred")) +
geom_text(aes(label=count)) +
theme(text=element_text(size =15), axis.text.x = element_text(angle = 90, size = 15), legend.position = "bottom")# Total richness
tar_load(faunaclass)
total_rich <- nrow(faunaclass %>%
distinct(., Species)) #725 species
# Total richness by sediment position
epi_rich <- nrow(faunaclass %>%
filter(Position == "Epifauna") %>%
distinct(., Species)) #341 epifaunal species
inf_rich <- nrow(faunaclass %>%
filter(Position == "Infauna") %>%
distinct(., Species)) #266 infaunal species
int_rich <- nrow(faunaclass %>%
filter(Position == "Interstice") %>%
distinct(., Species)) #118 interstitial speciesThere are 725 species in total, 341 being epifaunal, 266 being infaunal, and 118 interstitial.
tar_read(relab_phy)
tar_read(relab_ord) p.s.: it is normal to have NAs for some orders as higher classification for some of the species is ambiguous or unknown.
tar_read(relab_arth)
tar_read(relab_ann)
tar_read(relab_mol)
tar_read(correl)tar_read(densselR2) variables order R2 R2Cum AdjR2Cum F pvalue
1 Current_mean 6 0.17694646 0.1769465 0.1745677 74.385774 0.001
2 Gravel 2 0.04042285 0.2173693 0.2128323 17.819240 0.001
3 OM 3 0.03316824 0.2505375 0.2440015 15.224076 0.001
4 Fetch_max 8 0.05378450 0.3043220 0.2962092 26.518138 0.001
5 PC1_score 9 0.02568327 0.3300053 0.3202101 13.110071 0.001
6 T_mean 4 0.01314023 0.3431455 0.3315880 6.821629 0.012
tar_load(c(modfixdens12, modfdsel, modrdsel))
jtools::export_summs(
modfixdens12,
modfdsel,
modrdsel,
model.names = c("Density ~ Complexity", "Density ~ Complexity + Environment", "Density ~ Complexity + Environment + (1|Site)"),
statistics = c(
N = "nobs",
DF = "df.residual",
AIC = "AIC",
BIC = "BIC",
`Marginal R2` = "r.squared.fixed",
`Conditional R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "PC1_score",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "PC1_score:PC2_score",
"Mean current velocity" = "Current_mean",
"Depth",
"Exposure" = "Fetch_max",
"Gravel",
"Organic Matter %" = "OM",
"Mean Temperature" = "T_mean",
"Year"
),
scale = TRUE,
error_pos = "right",
bold_signif = .05
)| Density ~ Complexity | Density ~ Complexity + Environment | Density ~ Complexity + Environment + (1|Site) | ||||
| Rhodolith complexity (PC1) | 0.42 *** | (0.05) | 0.21 *** | (0.06) | 0.06 | (0.05) |
| Bed complexity (PC2) | 0.09 | (0.05) | -0.07 | (0.05) | 0.13 | (0.07) |
| PC1:PC2 | -0.31 *** | (0.05) | -0.14 ** | (0.05) | 0.04 | (0.04) |
| Mean current velocity | -0.27 *** | (0.05) | -0.50 | (0.24) | ||
| Depth | -0.29 *** | (0.06) | 0.08 | (0.14) | ||
| Exposure | 0.20 ** | (0.07) | 0.18 | (0.18) | ||
| Gravel | 0.24 *** | (0.04) | 0.05 | (0.05) | ||
| Organic Matter % | 0.20 *** | (0.05) | 0.14 ** | (0.05) | ||
| Mean Temperature | 0.11 ** | (0.04) | 0.06 | (0.03) | ||
| Year | 0.03 | (0.04) | 0.05 | (0.03) | ||
| N | 348 | 348 | 348 | |||
| DF | 344.00 | 337.00 | 335.00 | |||
| AIC | 889.58 | 767.68 | 687.22 | |||
| BIC | 908.84 | 813.91 | 737.30 | |||
| Marginal R2 | 0.24 | |||||
| Conditional R2 | 0.26 | 0.50 | 0.69 | |||
| Adjusted R2 | 0.25 | 0.48 | ||||
| F | 39.42 | 33.22 | ||||
| p-value | 0.00 | 0.00 | ||||
| All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_load(c(modfdann, modrdann, modfdart, modrdart, modfdmol, modrdmol, modfdepi, modrdepi, modfdinf, modrdinf, modfdint, modrdint))
jtools::export_summs(
modfdann,
modrdann,
modfdart,
modrdart,
modfdmol,
modrdmol,
model.names = c("Annelida LM","Annelida LMM", "Arthropoda LM", "Arthropoda LMM", "Mollusca LM", "Mollusca LMM"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`Marginal R2` = "r.squared.fixed",
`Conditional R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "PC1_score",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "PC1_score:PC2_score",
"Mean current velocity" = "Current_mean",
"Depth",
"Exposure" = "Fetch_max",
"Gravel",
"Organic Matter %" = "OM",
"Mean Temperature" = "T_mean",
"Year"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| Annelida LM | Annelida LMM | Arthropoda LM | Arthropoda LMM | Mollusca LM | Mollusca LMM | |||||||
| Rhodolith complexity (PC1) | 0.05 | (0.05) | 0.07 | (0.05) | 0.31 *** | (0.07) | 0.03 | (0.05) | 0.02 | (0.05) | -0.02 | (0.06) |
| Bed complexity (PC2) | 0.11 * | (0.05) | 0.11 | (0.07) | 0.00 | (0.06) | 0.10 | (0.06) | 0.31 *** | (0.05) | 0.18 ** | (0.07) |
| PC1:PC2 | -0.01 | (0.04) | 0.05 | (0.05) | -0.18 *** | (0.05) | 0.06 | (0.04) | -0.05 | (0.04) | -0.01 | (0.05) |
| Mean current velocity | -0.29 *** | (0.06) | -0.10 | (0.22) | ||||||||
| Depth | -0.31 *** | (0.05) | -0.21 * | (0.10) | -0.42 *** | (0.06) | -0.36 *** | (0.10) | ||||
| Exposure | -0.29 *** | (0.06) | -0.37 ** | (0.12) | 0.30 *** | (0.07) | 0.31 | (0.18) | ||||
| Gravel | 0.31 *** | (0.05) | 0.02 | (0.04) | 0.17 *** | (0.04) | 0.05 | (0.05) | ||||
| Organic Matter % | 0.06 | (0.05) | 0.22 *** | (0.06) | 0.10 * | (0.05) | 0.09 * | (0.05) | 0.12 * | (0.06) | ||
| Mean Temperature | 0.09 * | (0.04) | 0.07 * | (0.04) | 0.10 ** | (0.04) | 0.09 * | (0.04) | ||||
| Year | 0.15 *** | (0.04) | 0.16 *** | (0.03) | 0.06 | (0.04) | 0.07 * | (0.03) | -0.06 | (0.04) | -0.06 | (0.04) |
| N | 348 | 348 | 348 | 348 | 348 | 348 | ||||||
| Residual DF | 338.00 | 338.00 | 339.00 | 337.00 | 339.00 | 337.00 | ||||||
| Marginal R2 | 0.47 | 0.04 | 0.35 | |||||||||
| Conditional R2 | 0.59 | 0.58 | 0.37 | 0.75 | 0.54 | 0.47 | ||||||
| Adjusted R2 | 0.58 | 0.36 | 0.53 | |||||||||
| F | 53.91 | 24.94 | 49.88 | |||||||||
| p-value | 0.00 | 0.00 | 0.00 | |||||||||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||||||||
tar_read(anndenspc1)
tar_read(anndenspc2)
tar_read(artdenspc1)
tar_read(artdenspc2)
tar_read(moldenspc1)
tar_read(moldenspc2)Main phyla density as a function of HC
jtools::export_summs(
modfdepi,
modrdepi,
modfdinf,
modrdinf,
modfdint,
modrdint,
model.names = c("Epifauna LM", "Epifauna LMM", "Infauna LM", "Infauna LMM", "Interstitial fauna LM", "Interstitial fauna LMM"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`Marginal R2` = "r.squared.fixed",
`Conditional R2` = "r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "PC1_score",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "PC1_score:PC2_score",
"Mean current velocity" = "Current_mean",
"Depth",
"Exposure" = "Fetch_max",
"Gravel",
"Organic Matter %" = "OM",
"Mud",
"Mean Temperature" = "T_mean",
"Temperature Variation (sd)" = "T_sd",
"Year"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| Epifauna LM | Epifauna LMM | Infauna LM | Infauna LMM | Interstitial fauna LM | Interstitial fauna LMM | |||||||
| Rhodolith complexity (PC1) | 0.31 *** | (0.06) | 0.07 | (0.05) | 0.10 | (0.05) | 0.01 | (0.06) | 0.06 | (0.05) | 0.04 | (0.05) |
| Bed complexity (PC2) | 0.01 | (0.05) | 0.11 | (0.07) | 0.17 *** | (0.05) | 0.14 | (0.07) | 0.22 *** | (0.04) | 0.20 ** | (0.06) |
| PC1:PC2 | -0.09 | (0.05) | 0.07 | (0.04) | -0.11 * | (0.05) | -0.00 | (0.05) | -0.09 * | (0.04) | -0.03 | (0.05) |
| Mean current velocity | -0.28 *** | (0.06) | -0.23 | (0.21) | ||||||||
| Depth | -0.48 *** | (0.06) | -0.20 | (0.12) | ||||||||
| Exposure | 0.27 *** | (0.07) | 0.21 | (0.18) | -0.34 *** | (0.06) | -0.40 *** | (0.09) | ||||
| Gravel | 0.38 *** | (0.05) | 0.07 | (0.05) | -0.25 *** | (0.05) | -0.10 | (0.05) | 0.26 *** | (0.04) | 0.19 *** | (0.05) |
| Organic Matter % | 0.27 *** | (0.05) | 0.16 ** | (0.05) | 0.03 | (0.05) | 0.04 | (0.05) | ||||
| Mud | 0.33 *** | (0.05) | 0.25 *** | (0.06) | ||||||||
| Mean Temperature | 0.08 | (0.04) | 0.04 | (0.03) | 0.08 | (0.04) | 0.05 | (0.04) | ||||
| Temperature Variation (sd) | -0.01 | (0.04) | -0.01 | (0.06) | ||||||||
| Year | 0.00 | (0.04) | 0.02 | (0.03) | 0.09 * | (0.04) | 0.09 * | (0.04) | 0.14 *** | (0.03) | 0.14 *** | (0.03) |
| N | 348 | 348 | 348 | 348 | 348 | 348 | ||||||
| Residual DF | 338.00 | 336.00 | 339.00 | 337.00 | 339.00 | 337.00 | ||||||
| Marginal R2 | 0.09 | 0.11 | 0.56 | |||||||||
| Conditional R2 | 0.42 | 0.66 | 0.43 | 0.49 | 0.60 | 0.62 | ||||||
| F | 27.06 | 32.32 | 64.41 | |||||||||
| p-value | 0.00 | 0.00 | 0.00 | |||||||||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||||||||
tar_read(epidenspc1)
tar_read(epidenspc2)
tar_read(infdenspc1)
tar_read(infdenspc2)
tar_read(intdenspc1)
tar_read(intdenspc2)Density of macrofauna as a function of HC at each compartment
Figure 2
tar_read(psdens)
tar_read(psdtrait)Figure 2 (main text)
Check if the relationships between HC and species richness are the same for the total macrofauna as well as for the main phyla and the main compartments (different positions in sediment) separately.
Now we first forward select environmental variables and then add sites as a random factor to control for pseudo-replication and other site-related variables that were not taken into account.
Figure 3
tar_read(totrichpc1)
tar_read(figure3)Figure 3. OLS linear regression of Total Macrofauna Richness as a function of HC
tar_read(richselR2) variables order R2 R2Cum AdjR2Cum F pvalue
1 Depth 7 0.329480820 0.3294808 0.3275429 170.018051 0.001
2 Fetch_max 8 0.049649871 0.3791307 0.3755314 27.589068 0.001
3 T_mean 4 0.019017288 0.3981480 0.3928993 10.869694 0.002
4 Year 11 0.008659246 0.4068072 0.3998895 5.007009 0.023
5 PC2_score 10 0.008178651 0.4149859 0.4064330 4.781250 0.033
6 PC1_score 9 0.011702850 0.4266887 0.4166011 6.960742 0.005
7 OM 3 0.008783926 0.4354727 0.4238500 5.290328 0.017
tar_load(c(modfixrich12, modfrsel, modrrsel))
jtools::export_summs(
modfixrich12,
modfrsel,
modrrsel,
model.names = c("S ~ Complexity", "S ~ Complexity + Environment", "S ~ Complexity + Environment + (1|Site)"),
statistics = c(
N = "nobs",
DF = "df.residual",
AIC = "AIC",
BIC = "BIC",
`Marginal R2` = "r.squared.fixed",
`Conditional R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "PC1_score",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "PC1_score:PC2_score",
"Depth",
"Exposure" = "Fetch_max",
"Mean bottom temperature" = "T_mean",
"Organic Matter %" = "OM",
"Year"
),
scale = TRUE,
error_pos = "right",
bold_signif = .05
)| S ~ Complexity | S ~ Complexity + Environment | S ~ Complexity + Environment + (1|Site) | ||||
| Rhodolith complexity (PC1) | 8.38 *** | (0.91) | 2.89 * | (1.20) | 2.08 | (1.30) |
| Bed complexity (PC2) | 8.03 *** | (0.92) | 4.30 *** | (1.15) | 4.05 * | (1.58) |
| PC1:PC2 | -0.87 | (0.97) | 0.04 | (0.91) | 2.04 | (1.10) |
| Depth | -5.33 *** | (1.25) | -3.06 | (2.45) | ||
| Exposure | -2.75 | (1.41) | -1.54 | (2.83) | ||
| Mean bottom temperature | 2.05 * | (0.88) | 1.38 | (0.86) | ||
| Organic Matter % | 2.56 * | (1.12) | 3.72 ** | (1.23) | ||
| Year | 1.71 | (0.87) | 1.78 * | (0.83) | ||
| N | 348 | 348 | 348 | |||
| DF | 344.00 | 339.00 | 337.00 | |||
| AIC | 2956.83 | 2905.35 | 2872.96 | |||
| BIC | 2976.10 | 2943.88 | 2915.33 | |||
| Marginal R2 | 0.29 | |||||
| Conditional R2 | 0.33 | 0.44 | 0.43 | |||
| Adjusted R2 | 0.32 | 0.42 | ||||
| F | 55.56 | 32.69 | ||||
| p-value | 0.00 | 0.00 | ||||
| All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
You might notice that PC1 was not selected as it has a very small effect size and doesn’t contribute enough to the adjusted R2. I kept it in the final model anyways for comparison with the complexity-only model.
tar_load(c(modfixrich1, modfixrichpoly1, modfixrich2, modfixrichpoly2, modfixrich12, modfixrichpoly12))
jtools::export_summs(
modfixrich1,
modfixrichpoly1,
modfixrich2,
modfixrichpoly2,
modfixrich12,
modfixrichpoly12,
# model.names = c("Linear", "Quadratic"),
statistics = c(
N = "nobs",
DF = "df.residual",
AIC = "AIC",
BIC = "BIC",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
# coefs = c(
# "Rhodolith complexity (PC1)" = "PC1_score",
# "Bed complexity (PC2)" = "PC2_score",
# "PC1:PC2" = "PC1_score:PC2_score",
# "Depth",
# "Exposure" = "Fetch_max",
# "Mean bottom temperature" = "T_mean",
# "Organic Matter %" = "OM",
# "Year"
# ),
scale = TRUE,
error_pos = "right",
bold_signif = .05
)| Model 1 | Model 2 | Model 3 | Model 4 | Model 5 | Model 6 | |||||||
| (Intercept) | 68.07 *** | (1.00) | 68.07 *** | (1.00) | 68.07 *** | (1.00) | 68.07 *** | (1.00) | 68.07 *** | (0.90) | 66.67 *** | (0.96) |
| PC1_score | 8.22 *** | (1.00) | 8.38 *** | (0.91) | ||||||||
| `poly(PC1_score, 2)`1 | 8.21 *** | (1.00) | 6.96 *** | (1.02) | ||||||||
| `poly(PC1_score, 2)`2 | 1.37 | (1.00) | -1.25 | (1.06) | ||||||||
| PC2_score | 8.14 *** | (1.00) | 8.03 *** | (0.92) | ||||||||
| `poly(PC2_score, 2)`1 | 8.14 *** | (1.00) | 8.89 *** | (1.17) | ||||||||
| `poly(PC2_score, 2)`2 | 0.44 | (1.00) | 1.30 | (1.20) | ||||||||
| PC1_score:PC2_score | -0.87 | (0.97) | ||||||||||
| `poly(PC1_score, 2)`1:`poly(PC2_score, 2)`1 | 1.64 | (1.30) | ||||||||||
| `poly(PC1_score, 2)`2:`poly(PC2_score, 2)`1 | 4.02 ** | (1.25) | ||||||||||
| `poly(PC1_score, 2)`1:`poly(PC2_score, 2)`2 | -1.80 | (1.30) | ||||||||||
| `poly(PC1_score, 2)`2:`poly(PC2_score, 2)`2 | -3.58 ** | (1.19) | ||||||||||
| N | 348 | 348 | 348 | 348 | 348 | 348 | ||||||
| DF | 346.00 | 345.00 | 346.00 | 345.00 | 344.00 | 339.00 | ||||||
| AIC | 3028.54 | 3028.65 | 3029.73 | 3031.53 | 2956.83 | 2945.82 | ||||||
| BIC | 3040.09 | 3044.06 | 3041.28 | 3046.94 | 2976.10 | 2984.35 | ||||||
| Adjusted R2 | 0.16 | 0.16 | 0.16 | 0.16 | 0.32 | 0.35 | ||||||
| F | 67.23 | 34.64 | 65.82 | 32.93 | 55.56 | 24.45 | ||||||
| p-value | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | ||||||
| All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||||||||
Figure 4
tar_read(psrich)tar_load(c(modrrann, modrrart, modrrmol, modrrepi, modrrinf, modrrint))
jtools::export_summs(
modrrann,
modrrart,
modrrmol,
model.names = c("Annelida", "Arthropoda", "Mollusca"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`Marginal R2` = "r.squared.fixed",
`Conditional R2` = "r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "PC1_score",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "PC1_score:PC2_score"
),
scale = TRUE,
error_pos = "right",
bold_signif = .05
)| Annelida | Arthropoda | Mollusca | ||||
| Rhodolith complexity (PC1) | 1.50 * | (0.66) | 0.43 | (0.44) | -0.15 | (0.48) |
| Bed complexity (PC2) | 2.50 ** | (0.82) | 0.90 | (0.54) | 0.11 | (0.60) |
| PC1:PC2 | 1.50 ** | (0.56) | 0.59 | (0.36) | 0.06 | (0.40) |
| N | 348 | 348 | 348 | |||
| Residual DF | 342.00 | 342.00 | 342.00 | |||
| Marginal R2 | 0.11 | 0.02 | 0.00 | |||
| Conditional R2 | 0.38 | 0.52 | 0.51 | |||
| p-value | ||||||
| All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(annrichpc1)
tar_read(annrichpc2)
tar_read(artrichpc1)
tar_read(artrichpc2)
tar_read(molrichpc1)
tar_read(molrichpc2)Main phyla richness as a function of HC
jtools::export_summs(
modrrepi,
modrrinf,
modrrint,
model.names = c("Epifauna", "Infauna", "Interstitial fauna"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`Marginal R2` = "r.squared.fixed",
`Conditional R2` = "r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "PC1_score",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "PC1_score:PC2_score"
),
scale = TRUE,
error_pos = "right",
bold_signif = .05
)| Epifauna | Infauna | Interstitial fauna | ||||
| Rhodolith complexity (PC1) | 0.78 | (0.74) | 0.37 | (0.56) | 0.80 ** | (0.29) |
| Bed complexity (PC2) | 1.47 | (0.93) | 1.35 | (0.70) | 1.35 *** | (0.35) |
| PC1:PC2 | 1.40 * | (0.62) | 0.76 | (0.47) | 0.36 | (0.25) |
| N | 348 | 348 | 347 | |||
| Residual DF | 342.00 | 342.00 | 341.00 | |||
| Marginal R2 | 0.03 | 0.03 | 0.15 | |||
| Conditional R2 | 0.52 | 0.44 | 0.31 | |||
| p-value | ||||||
| All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(epirichpc1)
tar_read(epirichpc2)
tar_read(infrichpc1)
tar_read(infrichpc2)
tar_read(intrichpc1)
tar_read(intrichpc2)Species richness as a function of HC at each faunal compartment
As we can see, the marginal R2 for the models containing only a subset of the macrofauna are quite low. However, although the effects of HC are less evident in these cases, The pattern for most subsets is very similar to the general patter, with very slight positive effect of rhodolith complexity and a more important positive effect of bed complexity. Nevertheless Rhodolith complexity seems to play a more important role in driving annelid richness than the other main phyla.
tar_read(psN2)tar_read(psJ)tar_read(spfauna)Final PCA showing only the species with a godness of fit \(\geq\) .3
Figure 5
tar_read(faunapca)tar_read(faunaclust)After variable selection, we add sites as a factor to understand if we describe them well with the environmental variables we chose. We can look at the first two axes to have an idea of the relationships between the different explanatory variables.
tar_read(rdafsel) variables order R2 R2Cum AdjR2Cum F pvalue
1 Mud 2 0.087006599 0.0870066 0.08436789 32.973166 0.001
2 Depth 8 0.052115736 0.1391223 0.13413174 20.885579 0.001
3 Current_mean 7 0.052409189 0.1915315 0.18448093 22.299894 0.001
4 Branching_density 10 0.037531459 0.2290630 0.22007246 16.698239 0.001
5 Fetch_max 9 0.027419847 0.2564828 0.24561270 12.612470 0.001
6 Gravel 3 0.020457534 0.2769404 0.26421791 9.647916 0.001
7 T_sd 6 0.018777784 0.2957181 0.28121823 9.065187 0.001
8 Year 1 0.012098797 0.3078169 0.29148224 5.925444 0.001
9 L 11 0.010714212 0.3185312 0.30038554 5.314115 0.001
10 DR3 14 0.010718859 0.3292500 0.30934645 5.385398 0.001
11 T_mean 5 0.008804358 0.3380544 0.31638353 4.469044 0.001
12 OM 4 0.006918997 0.3449734 0.32150973 3.538580 0.001
13 Sphericity 13 0.005596140 0.3505695 0.32529227 2.878077 0.001
14 Total_Density 12 0.003994117 0.3545636 0.32742816 2.060685 0.001
tar_load(rdasite)
tar_read(triplotrda)RsquareAdj(rdasite)$adj.r.squared[1] 0.4279943
We check the model’s validity. Only 99 permutations were chosen for the validation by axis and by term here as they are really time consuming but the results are the same with 999 permutations. (Output is collapsed, click on code to see)
tar_read(aovsite)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + Total_Density + L + OM + Year + Site, data = envcomp)
## Df Variance F Pr(>F)
## Model 23 0.28507 12.289 0.001 ***
## Residual 324 0.32679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tar_read(aovaxis)
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 99
##
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + Total_Density + L + OM + Year + Site, data = envcomp)
## Df Variance F Pr(>F)
## RDA1 1 0.07620 75.5463 0.01 **
## RDA2 1 0.04629 45.8898 0.01 **
## RDA3 1 0.04244 42.0730 0.01 **
## RDA4 1 0.02841 28.1713 0.01 **
## RDA5 1 0.01876 18.5971 0.01 **
## RDA6 1 0.01339 13.2736 0.01 **
## RDA7 1 0.01225 12.1454 0.01 **
## RDA8 1 0.00957 9.4930 0.01 **
## RDA9 1 0.00781 7.7480 0.01 **
## RDA10 1 0.00710 7.0381 0.01 **
## RDA11 1 0.00557 5.5256 0.01 **
## RDA12 1 0.00307 3.0402 0.01 **
## RDA13 1 0.00239 2.3719 0.04 *
## RDA14 1 0.00202 2.0023 0.12
## RDA15 1 0.00176 1.7401 0.38
## RDA16 1 0.00154 1.5284 0.72
## RDA17 1 0.00133 1.3138 0.96
## RDA18 1 0.00113 1.1189 1.00
## RDA19 1 0.00103 1.0169 1.00
## RDA20 1 0.00088 0.8709 1.00
## RDA21 1 0.00077 0.7609 1.00
## RDA22 1 0.00073 0.7236 1.00
## RDA23 1 0.00065 0.6482 1.00
## Residual 324 0.32679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tar_read(aovterm)
## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 99
##
## Model: rda(formula = bcdens ~ Mud + Depth + Current_mean + Branching_density + Fetch_max + T_sd + Sphericity + T_mean + DR3 + Gravel + Total_Density + L + OM + Year + Site, data = envcomp)
## Df Variance F Pr(>F)
## Mud 1 0.05324 52.7815 0.01 **
## Depth 1 0.03189 31.6154 0.01 **
## Current_mean 1 0.03207 31.7934 0.01 **
## Branching_density 1 0.02296 22.7680 0.01 **
## Fetch_max 1 0.01678 16.6339 0.01 **
## T_sd 1 0.01202 11.9166 0.01 **
## Sphericity 1 0.00898 8.9014 0.01 **
## T_mean 1 0.00709 7.0263 0.01 **
## DR3 1 0.00610 6.0457 0.01 **
## Gravel 1 0.00858 8.5018 0.01 **
## Total_Density 1 0.00236 2.3400 0.01 **
## L 1 0.00491 4.8632 0.01 **
## OM 1 0.00426 4.2272 0.01 **
## Year 1 0.00573 5.6773 0.01 **
## Site 9 0.06813 7.5051 0.01 **
## Residual 324 0.32679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Figure 6
tar_read(upsetmain)tar_read(upsetsep)Figure 7
tar_read(bdplot)tar_read(replplot)tar_load(c(modbdtot, modrdftot, modrptot))
jtools::export_summs(
modbdtot,
modrptot,
modrdftot,
model.names = c("BDtotal", "Replacement", "Richness Difference"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
"Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
"PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| BDtotal | Replacement | Richness Difference | ||||
| Rhodolith complexity (PC1) | -0.51 ** | (0.16) | -0.42 * | (0.17) | -0.29 | (0.18) |
| Rhodolith complexity (PC1)^2 | -0.23 | (0.16) | -0.10 | (0.17) | -0.27 | (0.19) |
| Bed complexity (PC2) | -0.18 | (0.17) | -0.21 | (0.19) | -0.02 | (0.20) |
| PC1:PC2 | 0.28 | (0.19) | 0.11 | (0.21) | 0.32 | (0.22) |
| PC1ˆ2:PC2 | -0.07 | (0.20) | -0.27 | (0.22) | 0.24 | (0.24) |
| N | 30 | 30 | 30 | |||
| Residual DF | 24.00 | 24.00 | 24.00 | |||
| R2 | 0.43 | 0.32 | 0.20 | |||
| Adjusted R2 | 0.31 | 0.18 | 0.03 | |||
| F | 3.55 | 2.23 | 1.20 | |||
| p-value | 0.02 | 0.08 | 0.34 | |||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(bdtotpc1)
tar_read(bdtotpc2)
tar_read(rdtotpc1)
tar_read(rdtotpc2)
tar_read(rptotpc1)
tar_read(rptotpc2)tar_load(c(modbdann, modrdfann, modrpann))
jtools::export_summs(
modbdann,
modrpann,
modrdfann,
model.names = c("BDtotal", "Replacement", "Richness Difference"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
"Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
"PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| BDtotal | Replacement | Richness Difference | ||||
| Rhodolith complexity (PC1) | -0.61 *** | (0.14) | -0.64 *** | (0.15) | -0.16 | (0.19) |
| Rhodolith complexity (PC1)^2 | -0.27 | (0.15) | -0.11 | (0.15) | -0.30 | (0.19) |
| Bed complexity (PC2) | -0.08 | (0.16) | 0.01 | (0.16) | -0.14 | (0.20) |
| PC1:PC2 | 0.22 | (0.17) | 0.10 | (0.18) | 0.22 | (0.23) |
| PC1ˆ2:PC2 | -0.10 | (0.19) | -0.22 | (0.19) | 0.12 | (0.24) |
| N | 30 | 30 | 30 | |||
| Residual DF | 24.00 | 24.00 | 24.00 | |||
| R2 | 0.51 | 0.49 | 0.18 | |||
| Adjusted R2 | 0.41 | 0.38 | 0.01 | |||
| F | 4.98 | 4.62 | 1.03 | |||
| p-value | 0.00 | 0.00 | 0.42 | |||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(bdannpc1)
tar_read(bdannpc2)
tar_read(rdannpc1)
tar_read(rdannpc2)
tar_read(rpannpc1)
tar_read(rpannpc2)tar_load(c(modbdart, modrdfart, modrpart))
jtools::export_summs(
modbdart,
modrpart,
modrdfart,
model.names = c("BDtotal", "Replacement", "Richness Difference"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
"Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
"PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| BDtotal | Replacement | Richness Difference | ||||
| Rhodolith complexity (PC1) | -0.42 * | (0.17) | -0.29 | (0.16) | -0.38 * | (0.18) |
| Rhodolith complexity (PC1)^2 | -0.09 | (0.18) | -0.20 | (0.16) | 0.03 | (0.18) |
| Bed complexity (PC2) | -0.23 | (0.19) | -0.40 * | (0.18) | -0.02 | (0.20) |
| PC1:PC2 | 0.22 | (0.21) | -0.13 | (0.19) | 0.41 | (0.22) |
| PC1ˆ2:PC2 | -0.11 | (0.22) | -0.41 | (0.21) | 0.16 | (0.23) |
| N | 30 | 30 | 30 | |||
| Residual DF | 24.00 | 24.00 | 24.00 | |||
| R2 | 0.31 | 0.39 | 0.25 | |||
| Adjusted R2 | 0.16 | 0.27 | 0.09 | |||
| F | 2.14 | 3.11 | 1.57 | |||
| p-value | 0.10 | 0.03 | 0.21 | |||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(bdartpc1)
tar_read(bdartpc2)
tar_read(rdartpc1)
tar_read(rdartpc2)
tar_read(rpartpc1)
tar_read(rpartpc2)tar_load(c(modbdmol, modrdfmol, modrpmol))
jtools::export_summs(
modbdmol,
modrpmol,
modrdfmol,
model.names = c("BDtotal", "Replacement", "Richness Difference"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
"Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
"PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| BDtotal | Replacement | Richness Difference | ||||
| Rhodolith complexity (PC1) | -0.22 | (0.16) | -0.06 | (0.19) | -0.24 | (0.18) |
| Rhodolith complexity (PC1)^2 | -0.32 | (0.16) | -0.16 | (0.20) | -0.23 | (0.19) |
| Bed complexity (PC2) | -0.35 | (0.17) | -0.22 | (0.21) | -0.20 | (0.20) |
| PC1:PC2 | 0.30 | (0.19) | 0.19 | (0.23) | 0.15 | (0.22) |
| PC1ˆ2:PC2 | 0.15 | (0.21) | 0.02 | (0.25) | 0.20 | (0.24) |
| N | 30 | 30 | 30 | |||
| Residual DF | 24.00 | 24.00 | 24.00 | |||
| R2 | 0.40 | 0.13 | 0.19 | |||
| Adjusted R2 | 0.27 | -0.05 | 0.03 | |||
| F | 3.18 | 0.72 | 1.16 | |||
| p-value | 0.02 | 0.62 | 0.36 | |||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(bdmolpc1)
tar_read(bdmolpc2)
tar_read(rdmolpc1)
tar_read(rdmolpc2)
tar_read(rpmolpc1)
tar_read(rpmolpc2)tar_load(c(modbdepi, modrdfepi, modrpepi))
jtools::export_summs(
modbdepi,
modrpepi,
modrdfepi,
model.names = c("BDtotal", "Replacement", "Richness Difference"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
"Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
"PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| BDtotal | Replacement | Richness Difference | ||||
| Rhodolith complexity (PC1) | -0.42 * | (0.18) | -0.25 | (0.19) | -0.38 * | (0.18) |
| Rhodolith complexity (PC1)^2 | -0.16 | (0.18) | -0.06 | (0.19) | -0.19 | (0.19) |
| Bed complexity (PC2) | -0.04 | (0.19) | -0.07 | (0.21) | 0.02 | (0.20) |
| PC1:PC2 | 0.29 | (0.21) | 0.08 | (0.23) | 0.36 | (0.22) |
| PC1ˆ2:PC2 | -0.05 | (0.23) | -0.26 | (0.25) | 0.18 | (0.24) |
| N | 30 | 30 | 30 | |||
| Residual DF | 24.00 | 24.00 | 24.00 | |||
| R2 | 0.27 | 0.15 | 0.22 | |||
| Adjusted R2 | 0.12 | -0.03 | 0.06 | |||
| F | 1.82 | 0.82 | 1.38 | |||
| p-value | 0.15 | 0.55 | 0.27 | |||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(bdepipc1)
tar_read(bdepipc2)
tar_read(rdepipc1)
tar_read(rdepipc2)
tar_read(rpepipc1)
tar_read(rpepipc2)tar_load(c(modbdinf, modrdfinf, modrpinf))
jtools::export_summs(
modbdinf,
modrpinf,
modrdfinf,
model.names = c("BDtotal", "Replacement", "Richness Difference"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
"Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
"PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| BDtotal | Replacement | Richness Difference | ||||
| Rhodolith complexity (PC1) | -0.56 *** | (0.14) | -0.43 * | (0.17) | -0.28 | (0.18) |
| Rhodolith complexity (PC1)^2 | -0.35 * | (0.14) | -0.25 | (0.17) | -0.20 | (0.19) |
| Bed complexity (PC2) | -0.26 | (0.15) | -0.16 | (0.18) | -0.19 | (0.20) |
| PC1:PC2 | 0.18 | (0.16) | 0.09 | (0.20) | 0.17 | (0.22) |
| PC1ˆ2:PC2 | -0.09 | (0.18) | -0.18 | (0.22) | 0.13 | (0.24) |
| N | 30 | 30 | 30 | |||
| Residual DF | 24.00 | 24.00 | 24.00 | |||
| R2 | 0.57 | 0.33 | 0.19 | |||
| Adjusted R2 | 0.48 | 0.20 | 0.03 | |||
| F | 6.30 | 2.41 | 1.15 | |||
| p-value | 0.00 | 0.07 | 0.36 | |||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(bdinfpc1)
tar_read(bdinfpc2)
tar_read(rdinfpc1)
tar_read(rdinfpc2)
tar_read(rpinfpc1)
tar_read(rpinfpc2)tar_load(c(modbdint, modrdfint, modrpint))
jtools::export_summs(
modbdint,
modrpint,
modrdfint,
model.names = c("BDtotal", "Replacement", "Richness Difference"),
statistics = c(
N = "nobs",
`Residual DF` = "df.residual",
`R2` = "r.squared",
`Adjusted R2` = "adj.r.squared",
`F` = "statistic",
`p-value` = "p.value"
),
coefs = c(
"Rhodolith complexity (PC1)" = "`poly(PC1_score, 2)`1",
"Rhodolith complexity (PC1)^2" = "`poly(PC1_score, 2)`2",
"Bed complexity (PC2)" = "PC2_score",
"PC1:PC2" = "`poly(PC1_score, 2)`1:PC2_score",
"PC1ˆ2:PC2" = "`poly(PC1_score, 2)`2:PC2_score"
),
scale = TRUE,
transform.response = TRUE,
error_pos = "right",
bold_signif = .05
)| BDtotal | Replacement | Richness Difference | ||||
| Rhodolith complexity (PC1) | -0.33 | (0.16) | -0.46 * | (0.17) | 0.12 | (0.18) |
| Rhodolith complexity (PC1)^2 | -0.12 | (0.17) | 0.17 | (0.17) | -0.31 | (0.18) |
| Bed complexity (PC2) | -0.40 * | (0.18) | -0.26 | (0.18) | -0.16 | (0.20) |
| PC1:PC2 | 0.26 | (0.20) | 0.01 | (0.20) | 0.27 | (0.22) |
| PC1ˆ2:PC2 | -0.02 | (0.21) | -0.32 | (0.22) | 0.30 | (0.23) |
| N | 30 | 30 | 30 | |||
| Residual DF | 24.00 | 24.00 | 24.00 | |||
| R2 | 0.38 | 0.35 | 0.25 | |||
| Adjusted R2 | 0.25 | 0.21 | 0.09 | |||
| F | 2.92 | 2.55 | 1.57 | |||
| p-value | 0.03 | 0.05 | 0.21 | |||
| All continuous variables are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||||
tar_read(bdintpc1)
tar_read(bdintpc2)
tar_read(rdintpc1)
tar_read(rdintpc2)
tar_read(rpintpc1)
tar_read(rpintpc2)Figure 8
tar_read(bdtrait)
tar_read(bdphyl)Figure 8 (main text)